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An asymptotically exact many body theory for spin polarized interacting fermions in a one- 
dimensional harmonic atom trap is developed using the bosonization method and including back- 
' ward scattering. In contrast to the Luttinger model, backscattering in the trap generates one-particle 

potentials which must be diagonalized simultaneously with the two-body interactions. Inclusion of 
backscattering becomes necessary because backscattering is the dominant interaction process be- 
tween confined identical one-dimensional fermions. The bosonization method is applied to the cal- 
culation of one-particle matrix elements at zero temperature. A detailed discussion of the validity of 
the results from bosonization is given, including a comparison with direct numerical diagonalization 
in fermionic Hilbert space. A model for the interaction coefficients is developed along the lines of 
the Luttinger model with only one coupling constant K. With these results, particle densities, the 
Wigner function, and the central pair correlation function are calculated and displayed for large 
j^H fermion numbers. It is shown how interactions modify these quantities. The anomalous dimension 

I 1 of the pair correlation function in the center of the trap is also discussed and found to be in accord 

, with the Luttinger model. 
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PACS numbers: PACS numbers: 71.10.Pm, 05.30.Fk, 03.75.Fi 



I. INTRODUCTION 



The achievement of Bose-Einstein condensation in dilute ultracold gases jlj stimulated the theoretical interest in 
trapped fermionic many body systems 0.0, especially their superfluid prop erties pi I E B, Hi Recent experimental 
successes in obtaining degeneracy in three-dimensional Fermi vapors p^. ITOL ITU Il2l Il3j intensified the interest in 
confined Fermi gases. 

Using microtrap technology [T3. ITol fTijl, UtI ITsj j , it is conceivable to realize a neutral ultracold quantum gas of 
\ trapped quasi one-dimensional degenerate fermions. 
f* — ■ In many cases, identical spin polarized fermions experience only a weak residual interaction because s-wave scattering 
is forbidden. This makes the question of interactions in a one-component spin polarized fermion system somewhat 
academic. Possible exceptions are Feshbach resonance enhanced scattering between atoms ^| and electric dipole- 
^5 ' dipole interactions in the case of polar molecules [2(j . 

The confinement of a trapped ultracold gas can be realized by a harmonic potential. We have developed an 
asymptotically exact theory of interacting one-dimensional fermions confined to a harmonic trap. It is based on the 
bosonization method known from the theory of Luttinger liquids (cf. I2TI l2p . I23I |24|] ) and exploits the linearity of the 
energy spectrum of free oscillator states. The method was presented in |25j for the one-component gas with forward 
scattering and extended in |2a | to the case of two components. This model must be seen as a soft boundary alternative 
I i to the well studied case of one-dimensional interacting fermions confined by hard walls [23, |2l| 123, bMJ, |3l|, |32, t^j ■ 

Meanwhile, an investigation of the interaction matrix elements for one-dimensional identical fermions in the har- 
monic trap revealed (34[ that backscattering dominates forward scattering unless the pair potential is long ranged. 
Unlike in the case of the Luttinger model, backscattering cannot be taken into account by merely renormalizing the 
forward coupling constants. This is due to the one-branch structure of the present model, which generates one-particle 
potentials from the backward scattering when the latter is brought into the form of an effective forward scattering. 
These one-particle potentials must be diagonalized simultaneously with the two particle interactions from forward 
scattering. We solve this problem by supplementing the squeezing transformation with an appropriate displacement 
transformation. 

The paper is organized as follows: Section II describes the model and classifies the scattering processes for fermions 
in the one-dimensional harmonic trap. Section III gives the solution of the backscattering problem for the one-particle 
matrix elements using the bosonization method. Section IV discusses the validity of the bosonization scheme and 
compares the results with those from direct numerical diagonalization in fermionic Hilbert space. Sec. V presents 
results for several quantities of interest: Occupation probabilities, particle and momentum densities, and the central 
pair correlation function, which are all derivable from the one-particle matrix elements. We employ a model of 
the interaction coefficients developed in analogy with the Luttinger model. It is characterized by just one coupling 
constant K and by a small decay parameter r<l. This approach is described in the Appendix. 
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II. THEORY 
A. Description of the model 

We consider spin polarized fermionic atoms interacting via the pair interaction operator 

V = \ V ( m 'P^ ( l' n ) (4A)(cj,Cn)- (1) 

mnpq 

The fermions are confined to a highly anisotropic axially symmetric harmonic trap. The trap potential is 

V(x, y, z) = ^m A u 2 z 2 + \m A ul {x 2 + y 2 ) = V z (z) + V p (p). (2) 

The atom mass is denoted uia and z is the one-dimensional coordinate in the elongated axial direction of the trap. 
The trap frequencies are uj£ and lo±_ 3> loi- The quasi one-dimensional Fermi gas is characterized by N (< u>±/u£) 
identical fermions filling the first N one-particle levels 

hu n = hoJt(n + l/2), n = 0,l,... (3) 

of the one-dimensional harmonic potential V z , while the transverse part of each unperturbed wave function remains 
the transverse ground state ip±o(x)ip±o(y)- 

The unperturbed one-dimensional Fermi energy is 

e F = Huj e (N - 1) + -nut = hu; e (N - -). (4) 

The Fermi energy is only slightly smaller than the excitation energy hui± of the first excited transverse state in the 
case of N = N max = largest integer in cu±/ui£. We have in mind a situation when the filling factor F = N/N max is 
small enough so that the Fermi level is well below this excitation energy. The above assumption about the transverse 
wave function is then justified. Furthermore, there still exists a macroscopic number of possible excitations of the 
one-dimensional Fermi sea that do not violate this condition. 

The operators in & m and c q in Eq. Q denote fermion creation and destruction operators, respectively. They obey 
the fermionic algebra CmC^ + c^dm = <5 mj „. This ensures that each oscillator state with single particle wave function 

U^) = ^^e-^H n (az) (5) 

and energy e n — huj n is at most singly occupied. The intrinsic length scale of the system is the oscillator length 
I = or 1 where a is defined by a 2 — mAOJe/h. H n denotes a Hermite polynomial. The non- interacting Hamiltonian 
. u) n c^c n has linear dispersion and incorporates exactly the harmonic trapping potential. 
The interaction matrix elements Vim, p; q, n) in Eq. Q are calculated from effective one-dimensional pair potentials 
using the harmonic oscillator states JSJ. Thus each individual interaction matrix element contains information about 
the harmonic trap. 

In a degenerate Fermi system, the most relevant states are those near the Fermi energy, i.e. msspwqwnwA^. 
This limited number of interaction matrix elements can be further reduced by a classification scheme based on 
approximate momentum conservation in the center of the trap and near the Fermi energy. This is described in the 
following subsection. 

B. Classification of coupling coefficients 

We are interested in the modification of the zero temperature Fermi sea due to the interaction Q ■ This problem 
is different from that of scattering between excitations above the Fermi sea, which, in the present one-dimensional 
context, are probably not quasi particles. 
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First of all, we must classify the coupling coefficients in Eq. QJ. There are 24 coupling coefficients V(m,p; q, n) for 
any set of four distinct integers. Due to the symmetries (to <-> q), (n «-> p), and ([to, q] [n,p]) only three coupling 
coefficients remain. It is shown in |34| that among them, those with index combinations n = m + p — q, n = p + q — to, 
and n = m + q — p dominate, provided the fermion number N is large. It is here, where the fermionic nature of the 
quantum gas becomes relevant. 

We can thus write 

V(m,p;q,n) ->V a 6 

m—q,n—p 

+ V b 5 v + v c s (6) 

Other weakly inhomogcncous trapping potentials also give dominant sets of V a ,Vb and V c coupling constants but 
do not lead to a linear dispersion required for strict bosonization. Note that the coupling constants still depend on 
indices, e.g. V c — V(m,p; q,n = m + q — p). In the present model, this is further simplified to V c — V c (\p — q\). 

Qualitatively Eq. © can be understood as follows: The single particle states ip n well inside the trap are super- 
positions of plane wave states exp(ik n z) with k n = ±a \/2n + 1. For N 1, the relevant states are near the Fermi 
energy and thus \k n \ « Uf = a^/2N — 1. Here, kp denotes the Fermi wave number. 

According to Eq. |JT|1. incoming states {n, q} are transformed into {p, to} in the collision process and the momenta 
of these states are (approximately, because of the weakly inhomogeneous trapping potential) conserved. Denoting a 
state with k n « — kp by (—n), three distinct collision processes can be discriminated: 

{n, q} -> {p, m}, {n, (-q)} -> {p, (-to)}, {n, (-q)} {(-p),m}. (7) 

Processes with strict momentum conservation dominate: This explains the Kronecker symbols in the approximate 
relation JBJ). Momentum transfer in the first two cases is small, thus describing forward scattering. In the last case, 
the momentum transfer is about 2kp, which corresponds to backward scattering. The first two cases were considered 
in jlSE^I- The last one requires an extension of the bosonization method, which is the aim of the present paper. 

The couplings V a , Vb, and V c are the analogues of the Luttinger model couplings 54, 52, and g\, respectively 
[2ll I23I. Izil] . In contrast to the Luttinger case, it was found in j^H that in a gas of identical fermions confined to the 
harmonic trap, forward scattering is almost completely suppressed so that backscattering is the dominant interaction 
process, unless the pair interaction potential is of long range. This is essentially a consequence of the Fermi algebra. 
In the following, we will ignore V a and Vb completely though a more general treatment is possible. 

In restricting the full set of interacting matrix elements to a set of solvable interactions we define a simplified model 
which cannot fully represent the initial problem. 

For a number of properties such as anomalous dimensions, we can, however, expect universality in the sense of the 
Luttinger liquid phenomenology [21,23]. This is confirmed by our result for the one-particle correlation function (Sec. 
V B) which shows Luttinger liquid behavior in the center of the trap. 

III. BACKSCATTERING AND BOSONIZATION 

The treatment of forward scattering is described in [25L |2^j . We give here only the extensions necessary for the 
inclusion of the backscattering interaction coefficients 

V(m,p;q,n) = V c (\q ~ p\) 8 m+q ^ +p . (8) 
Substituting these interaction coefficients into Eq. (QJ and reordering operators gives 

V c = -\ E V M) + \ E V c (\v\)(cl +V c m _ v ). (9) 

We omitted terms which are proportional to the fermion number operator by setting V c (0) = 0. The second term is 
a one-particle operator, which appears due to the backscattering. The perturbation Eq. is exactly solvable. Thus 
there is no renormalization group flow of the coupling V c . 

The essential requirement of the bosonization method is the possibility to express all operators entirely in terms of 
density fluctuation operators. It is still met in the present case: Introducing the density fluctuation operators 
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Pi?) = ^2 d l+p 5 9> 



(10) 



or, more conveniently, canonical boson operators related to them by: 



Pip) 



it is found that bosonic commutation relations 



\p\ d\ p \, p < 0, 
VP %, p > 0, 



(11) 



are satisfied after introducing the anomalous vacuum (cf. e.g., |24Ll35)l. 
The bosonized form of the backscattering operator is 



(12) 



m>0 



V c (m) + + ! } vlmVMnn 



2 ^ 

m>0 



l 2m 



(13) 



It is seen that the two- par ticle interaction due to backward scattering is of the same form as the forward scattering 
operator Vb studied in |25j except for a sign change, V& — * — V c . This is in complete analogy to the Luttinger case. 
However, the remaining one-particle operator produces non-trivial changes. 
In order to diagonalize the total Hamiltonian 



H = huji ^2 m d\nd m + V c , 



(14) 



m>0 



we perform two canonical transformations: 



The first one 



S. 



So. 



(15) 



Si = exp 



cy y t Cm(/m frn ) 



2 

. m>0 



(16) 



is a kind of squeezing transformation and was used in |25j to diagonalize the two-particle interactions. It gives 



{sl/mSij = f m cosh Cm - fin sinh C„ 



(17) 



The second one 



S 2 = exp 



^ ] Vmifm fm) 



rn>Q 



(18) 



is a displacement in order to get rid of the terms linear in the d and <v operators. The total result is 



dm = f m COsh Cm - f m Sinh Cm + Vm exp(^Cm)- 

We find two diagonalization conditions: The standard one 



(19) 
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tanh 2Q n 



V c (m) 



and another one due to backscattering 



-V c (m/2) exp(~C m )/(2^me m ), m = 2n, 

0, to — 2n — 1 . 



The final form of the Hamiltonian is 



H = H + V c = m e m P m j m + const. . 



m>0 



The renormalized oscillator frequencies are 



and 



defines the dimcnsionlcss coupling constants 



exp(-Cm) = VK„ 



(20) 



(21) 



(22) 



(23) 



(24) 



I huj( + V c (to) 
hut - V c (m) 



(25) 



One-particle matrix elements 

We a pply the theory to the evaluation of the one-particle matrix elements (c\ n c n ). To this order, we follow the 
steps in |25j and introduce the bosonic field 



n— 1 * 



(26) 



which allows to express a particle number conserving bilinear product of auxiliary Fermi fields studied by Schonhammer 
and Meden IH 



l=-c 



(27) 



tpl(u)ip a (v) = G N (u- v) expj-z ^ (u) — P (v) |exp|-i cj>(u) - <j>{v) } 
The quantity Gat (it) in Eq. (|28|l is a distribution valued Fermi sum defined by 



(28) 



N-l 



G N {u) = e ~ a(u 



i=-c 



In order to evaluate expectation values of exponentials containing 0-operators, one must 



(29) 
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a) express the d-operators in terms of the free / and /^-operators, 

b) apply the bosonic Wick theorem. 

The Wick theorem refers to homogeneous linear combinations of / and /^-operators. Due to backscattering, the 
operator (f> contains a c-number part <p c , which must be treated separately: 



with 



<p c (u) = -iC(-u), (30) 



C{u) = ^2 6m exp(— 2imu). (31) 



The quantities are given by 



&m = mm \ = -V c (m) 2m . (32) 

V 2m 4m e 2m 

Following the steps in [2^ , the zero temperature expectation value of Eq. (|28[) becomes 

(4>t(u)i> a (v)) = G N (u - v) exp[-Wi(u,v)-W 2 (u,v)}, (33) 

with W\ given by 

^ 2 

W\(u, v) = > — fo m — a m cos m(u + v)]{l — cos m(u — v) \ . (34) 



m 

m>0 



The interaction parameters are 



a m = - sinh2C m , J m = sinh Cm- (35) 
The contribution from the one-particle operator is: 

oo 

W 2 (u,v) = C(-u) - C{u) + C(v) - C(-v) = Ai V] £ 2m sinm(u - v) cosm(u + v). (36) 

ra— 1 

In order to obtain the matrix elements, a number of calculational steps are done, starting with the coordinate 
transformation (u + v)/2 = t, u — v = s. Using Wi(u, v) — Wi(s,t) and Wi(s ± 2n,t) = Wi(s,t ± n) — Wi(s,t), 
(i — 1,2) and further symmetries, one obtains 



; /" r ds r - , , - , ( exo\is(n - N + 1)1 , 

(ci_ p c n+p ) = / — cospt ^ — exp[-^( S ,t/2) - W 2 (s,t/2)] ( i _ exp(j a - e ) ^ " |: ' 57) 



The distribution in curly brackets can be written as 

"cos(n — N + h)s 



Jexp[is(n- 7V + 1)]1 J sin(n - TV + ±)s 
\ 1 - exp(is - e) J \ 2sin(s/2) 



+ i 



2sin(s/2) 



+ tt<5 2t ( S ). (38) 



By defining 



W 2 (s,t/2)=if(s,t), 



(39) 
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the final result becomes 




/. 



— TV 




(40) 



x 



— TV 



•TT 




sm[(n- N + l/2)s- f(s,t)} 
2sin(s/2) 



} 



exp[-Wi(s,t/2)]. 



Besides fixing the transformation parameters £ m in Eq. (1201) , the effect of backscattering appears in the argument 
of the sine via the function f(s,t), which is explicitly given by 



Backscattering thus destroys the specific form of particle-hole symmetry found in 25] for forward scattering. 



We can state that our solution Eq. (|40|) is the exact result for N one-dimensional fermions confined to a harmonic 
trap, interacting via Eq. lO, and immersed in an anomalous vacuum of fermions filling all negative energy states. 
Thus we have to assess the role of the anomalous vacuum for the present finite size system. Intuitively, it is clear that 
its role decreases with increasing Fermi energy ej? oc iV, the energy region where interactions are most relevant. 

We must, however, also consider the strength of the interaction. The dominant dimcnsionlcss coupling constant 
in Eq. I|25|l is K\ = K, since K m decreases with increasing m. It is then seen that K varies from zero to infinity 
when the physical coupling coefficient V c (l) varies from —Tiuji to hajg. Values of V c (l) outside this ran ge a r e p hysicall y 
inaccessible, as are coupling constants g = g 2 = g^ < —irTivp in the corresponding Luttinger model |2lT 1231 Y24 . l35 | . 
Note that at the extreme values V c (l) = ±TiuJt the renormalized excitation energy e\ according Eq. i|23[l vanishes. 

In the Luttinger model, K — > oo or g — > —nhvp corresponds to the strongest physically allowed attraction. The 
compressibility vanishes and a phase separation occurs (cf. the discussion in [23j). 

Our numerical results for the particle density in the present model show an increasing extension beyond the classical 
turning points, i.e., the density progressively leaks out of the trap for increasing K 3> 1. The interpretation of this 
effect is difficult due to the presence of the anomalous vacuum and the finiteness of TV: It is conceivable that an 
increasing interaction pulls more and more fermions out of the anomalous vacuum as indicated by studies of the 
fermion sum rule below. 

Nevertheless, we can make the following statement about the validity of our bosonization scheme: For any fixed 
|y c |/(?iwf) < 1, i.e., < K < oo, the error due to the anomalous vacuum can be made as small as wanted by increasing 
the physical particle number N. We call this asymptotically exact. 

Though we do not have an analytic expression for the error at present, we know from our study of the fermion sum 
rule and the results presented below that the error decays fast with increasing N, probably in an exponential way. 

The precise nature of the singular points K — and K ~ oo in the present finite system requires, however, further 
investigation. In assessing its relevance for interacting confined fermions, one must also consider the dependence of 
V c on the particle number. This in turn depends on the specific form of the interaction potential (cf. 34]). It thus 
seems that there is no simple and general limit N — > oo accompanied by a proper scaling of loi and the interaction 
parameters ("thermodynamic limit") for the interacting system. 

We expect, however, that the region near the center of the trap (|z| <C hp oc y/~N) acquires properties of a 
homogeneous Luttinger liquid. This is demonstrated for the central pair correlation function studied in the following 
Section and for a reasonably large particle number N = 10 3 . 

We also note that numerical investigation of the interaction coefficients in [34| show that increasing N does not 
alter the dominance of V c over the forward scattering coefficients in a gas of identical fermions. 




(41) 



m— 1 



IV. VALIDITY OF THE THEORY 



A. 



Numerical Method 



For the purpose of the present investigation, it is sufficient to use a simplified model ("toy" model IM1 in |2j|) by 
retaining only the terms with v = ±1 in Eq. The relevant parameter then are: 
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ai = , 71 = ^ — . (42) 

The function W%(s, t/2) in Eq. (|40|l becomes 2 [71 — aiCOs(t)][l — cos(s)] and the backscattering expression l|41|) 
simplifies to /(s,t) = — V c (l) sin(s) cos(t)/hLUg. The one-particle matrix elements at zero temperature are then given 
by 

t .... N- 1 ^ f d * ^ r ds { sin[(M + 1/2 -N) S -f(s,t)] 



[ c k _ p c M+p ) = -S p , -J^ -cosipt) — — — I (43) 

x exp{— 2 [71 — ot\ cos(t)][l — cos(s)]} . 

This is the prediction of the bosonization method for identical fermions in the one-dimensional harmonic trap with 
dominant backscattering. 

In |25j, it was already pointed out that first order perturbation theory in fermionic Hilbert space reproduces the 
results using the equation corresponding to Eq. (|4*3)l for forward scattering (f=0) when it is expanded to first order 
in o,\ and /3i, and the same applies to the present case: Both approaches lead to the weak coupling result 



(4- p cm+p) = <5 P ,o &{N - M - 1/2) - Sm,n-i {S p .i + 5,,-x) + 0^ 



hut 



(44) 



Finally, the predictions of Eq. I|43(l will be checked against numerical results obtained by direct diagonalization in 
the fermionic A^-particle Hilbert space and for strong coupling. 

The fermionic iV-particle Hilbert space is spanned by the unperturbed A^-particle states 

\{m})^^cl l cl 2 ---cl N \vac). (45) 

Here, {to} denotes a sequence of occupation numbers m n — 0, 1 for the single particle states tp n (z). For example, the 
three-particle state with energy (11/2) hu>? and excitation energy Tuog is 

|1,1,0, 1,0,0,... ) (0) = clc\cl\vac). (46) 

In order to simplify the notation for the unperturbed A^-particle states, we classify them according to their excitation 
energies AE(n) — nhu>£. Degeneracy is taken care of by ordering the states according to the lowest unoccupied single 
particle state occurring in them. They are then numbered consecutively \m)(°\m = 0,1,2,.... The unperturbed 
ground state is thus 

|0)<°> = |1,1,. ..,1,1,0,0,...), (47) 

and subsequent excited states are 

|1) (0) = |1,1,.. .,1,0,1,0,0,...), (48) 

|2)(°) = |1,1,... ,0,1,1,0,0,...), 

|3)<°> - |1,1,.. .,1,0,0,1, 0,0,...), etc.. 
The perturbed ground state to first order used to obtain Eq. (|4^|) is 

| )M = |0)(°) + (^r) |2)(°). (49) 
V 2 ^ / 
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This result was also used to check the numerical procedure. 

The actual eigenstates of the interacting Hamiltonian are denoted by |s). They are expanded according to 

|s>=5> S m|m)W (50 ) 

m 

These expansion coefficients are the central quantities in the numerical procedure: In terms of the expansion coeffi- 
cients, the occupation probabilities for oscillator states ipM with respect to the ground state |0) are 

P(M) = comC* 0n ^(nlcljCMlm)^. (51) 

van 

Similarly, the particle density, which is the expectation value of the density operator ip^ (z)ip(z) with 'ip(z) = 
^2 n i>n(z)c n , becomes 

OC 

n(z) = (0|^( Z )^( Z )| ) = Y, Mz)Mz)Y, c ° m $n (0 Hn\clc q \m)W. (52) 
The first computational step is the evaluation of the interaction matrix elements 

Vij = {0 Wc\j) {0) . (53) 

In detail, the free TV-particle eigenstates Eq. <|45ll and the expectation values in Eas. (|51(l . (|52|l . and (|53|l with respect to 
these states are calculated via an algorithm which implements the fermion-algebra. The number of eigenstates which 
are taken into account grows strongly with the maximal excitation energy AE(n max ) due to increasing degeneracy. 
In subsequent steps, the diagonalization of the matrix 

//,, //■„., • \ . r (54) 

and the computation of the expansion coefficients c sm follows j^E ■ Eigenvalues and eigenvectors of the matrix 
(|54|l are computed using a QL algorithm. The normalized eigenvector of the ground state determines the required 
expansion coefficients co m according to (°)(m|0). 

B. Occupation Probabilities 

The occupation probabilities of oscillator states 4>m in the interacting ground state are 

< P(M) = (0|4<W|0) < i. (55) 



Due to the backward scattering, they do not show the symmetry found in 25, 26] for forward scattering. 
We will also discuss the sum rule 



S(N, ai )=Y(^m)=N (56) 

m=0 

for the fermion number. This sum rule is only asymptotically fulfilled in the bosonization method: The particles in 
the anomalous vacuum couple to the physical particles, giving an effective particle number S(N,a\) > N .38]. The 
effect is most pronounced when the particle number N is smaller than the coupling strength |ai|. The coupling to 
the anomalous vacuum also leads to the surprising feature that interaction effects occur for just one physical particle 
(N = l). 
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2 4 6 8 10 12 14 16 18 20 



n 

Fig. 1: Occupation probabilities P(n) of oscillator states ip n in the interacting system for repulsive backscattering 
strength a\ — 1 at zero temperature and particle numbers N = 5 and N = 14. Crosses are results of numerical 
diagonalization, predictions of the bosonization method are given as squares (N = 5) and circles (N = 14), 
respectively. No significant deviations between the two approaches are seen for particle number N = 14. For 
N = 5, the presence of the anomalous vacuum in the bosonization method causes visible deviations. 

Figure 1 shows a comparison of occupation probabilities calculated using the bosonization method and numerical 
diagonalization. For a coupling constant a± = 1, which corresponds to V c {\) — — 0.894 tiwg, no significant deviations 
are found for particle numbers N = 14 and larger. However, for N = 5, the bosonization method is less accurate 
because the anomalous vacuum is present. Though the excess AN (5, 1) = 5(5, 1) — N amounts to only 5.4 • 10~ 4 , 
individual occupation probabilities are less accurate. 



C. Particle density 



The density of the trapped particles can be written as 



oo oo M 

z) (c]y _„CM-\-p)- (57) 

M=0 Af=lp=l 

The particle density shows Friedel oscillations |40l Elj . The amplitudes of the Friedel oscillations are modified by 
the interaction as described in |25J, 143J . 

Figure 2 displays particle densities for the coupling value a% = 10, corresponding to V c (l) = — 0.998 fiwe, which 
is a rather large repulsive interaction. The results of numerical diagonalization and of the bosonization method are 
compared for various particle numbers N. The deviations found are due to the presence of the anomalous vacuum 
in the bosonization method which interacts with the physical particles. The differences are more pronounced for 
larger coupling strength or smaller particle number N: Larger particle numbers lead to larger Fermi energies cf- 
This decouples the physically relevant energy region e w €f from the anomalous vacuum. Larger coupling strengths 
counteract the decoupling. 
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2 4 6 8 10 

Z/I 

Fig. 2: Dimensionless particle density n(z)l in units of the inverse of the oscillator length I versus dimensionless 
distance z/l from the center of the one-dimensional harmonic trap for several particle numbers N of a Fermi 
gas with dominant backward scattering at zero temperature. Thick curves are from numerical diagonalization 
and dotted curves show the Fricdcl oscillations calculated by the bosonization method. Interaction strength is 
a\ — 10. Deviations for particle numbers less than N = 30 are due to the anomalous vacuum, which couples 
strongly to the real particles. 

We can summarize: The results of the bosonization method for the one-particle matrix elements of a harmonically 
trapped one-dimensional Fermi gas with backscattering when compared with the exact numerical diagonalization 
show noticeable deviations from the exact results when the physical number N of fermions is small or the interaction 
strength is large. This effect is not restricted to backscattering. Forward scattering would produce similar results. For 
sufficiently large N, the bosonization method produces, however, fully acceptable results. It is then far more effective 
than any numerical approach. 



V. RESULTS FOR LARGE N 

The functions Wi(s, t) and f(s, t) in the exact expression (|40|) for the matrix elements contain a bewildering number 
of interaction coefficients, which all depend on the index m. In order to proceed with the evaluation, we describe in 
the Appendix a model for the interaction coefficients, which allows to do all summations and express them by just 
two parameters: The main coupling constant K and a small number r< 1 specifying the exponential decay of all 
interaction coefficients. Following [2(J, we adopt r = l/y/~N so that only the coupling constant K > remains. K > 1 
corresponds to attractive interactions and < K < 1 to repulsion between the fermions. An analogous procedure is 
used in the theory of the Luttinger model PII2II23 I l24 | . 

The result for the one-particle matrix elements becomes 



(c ] c + ) - ^ 



w dt 



cos(pt) 



w ds fsin[(n- N + l/2)s - f(s,t)} 



„ 2ir [1 + Z a - cos(t)] Q o 2ir 



[1 + Z~- cos(s)] 



2sin(s/2) 

\ {[1 + Z a - cos(t - a)] [1 + Z a - cos{t + s)]} ao/2 , 



(58) 
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with constants Z 7 and Z a according to 



cosh(r 7 ) - 1 -> r 2 /2, Z a = cosh(r Q /2) - 1 -> r 2 



(59) 



The coupling dependent exponents in Eq. I|58[) are 



(60) 



The sign of ao is negative when the backscattering is attractive. 

In addition, we take advantage of the fact that backscattering is dominant in a harmonically confined gas of identical 
fcrmions: The function / is then determined by the backscattering contribution C{u) in Eq. (|A13J) : 



f(s, t) = 23[C(</2 - s/2) - C(t/2 + s/2)] = - (1 - K 2 ) arctan 



2 q (cos £ — q cos s) sin s 



1 — q 2 — 2q (cos £ — q cos s) cos s 



(61) 



with 



q = exp(-r/2) -> 1 - r/2 



(62) 



A. Particle and momentum densities 



The particle density n{z) is evaluated using Eq. I|57|l and the matrix elements Eq. I|58|l . For large N and moderate 
coupling (K is of order unity), the particle density well inside the trap is dominated by the large number of nearly 
filled states inside the Fermi sea. However, the states above the Fermi level Np = N— 1 make significant contributions 
to n(z) by modifying the Friedel oscillations |4pJ: This is demonstrated in Fig. 3 using the full expression H58|) . 



14.30- 








43.0 43.5 44.0 44.5 45.0 45.5 



z/l 



Fig. 3: Dimensionless particle density n(z)l in units of the inverse of the oscillator length I versus dimensionless 
distance z/l near the classical boundary at Lp = 1\J2N — 1 of the one-dimensional harmonic trap for N = 1000 
interacting spinless fermions at zero temperature. Dotted curve shows unperturbed Friedel oscillations. Strongly 
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oscillating curve refers to a repulsive interaction with K = 1/3 while the dashed curve is for an attractive 
interaction with K = 3. The inset displays the density near the center of the trap, where Friedel oscillation 
have the usual period n/kp = wl/\/2N — 1. 

It is seen that Friedel oscillations are enhanced in the repulsive case while an attractive interaction suppresses 
them. Near the classical boundary \z\ w Lp = I y/2N — 1, the oscillation period is larger than the standard value 
7r jkp well inside the trap and obeys the relations found in |4*H for the non- interacting case. Furthermore, it is seen 
that the perturbed density extends considerably into the classically forbidden region. The effect becomes stronger 
with increasing coupling strength [K 1 or <C 1). 

The momentum density distribution is given by 



P( k ) = ' 2 EE (-l) p ^ n - P (kl 2 )^ n+P (kl 2 ) (ct_ p c n+p ). (63) 

n=0 p=—n 

It is seen that p(k) is identical in shape to the density distribution n{z) provided the diagonal approximation for 
the one-particle matrix elements is reasonable. This is the case for moderate attractive coupling. 

We know [2j| that the Friedel oscillations in n(z) and p(k) behave oppositely: Attractive interactions increase them 
in the momentum density and decrease them in the particle density and vice versa for repulsion. A corresponding 
effect is seen in the Wigner function discussed below. 

B. Wigner function 

In terms of the local creation and annihilation operators t/>T (z) and ip(z), the static Wigner function of the many- 
fermion system is given (cf. |42j1 by 

/oo 
dCe-^{^{z-Cl2)^{z + Cl2)). (64) 
-oo 

The Wigner function of the one-dimensional Fermi gas with two components and forward scattering between the two 
components was studied in |43j. 

Transforming to the oscillator representation gives 



W(z,k)= ^ (4»Cn) fmn(z,k), (65) 
m,n— 

with expansion coefficients f mn according to 

dCe- ik < 7p m {z - (/2)i> n (z + C/2) = f nm {z,-k), (66) 

-oo 

These are explicitly given in terms of generalized Laguerre polynomials (n > m, cf. e.g., |44| ) by 

f mn (z, k) = 2(-l) m (2"- m m!/n!) 1 / 2 (z/l - ikl) n - ,n exp {-z 2 /l 2 - k 2 l 2 ) L^- m > (2z 2 /I 2 + 2k 2 l 2 ) . (67) 
It is noted that the one-particle matrix elements can be completely reconstructed from the Wigner function: 

1 f°° 

{& m c n ) = — / dzdk f nm (z,k)W(z,k). (68) 

The Wigner function is thus equivalent to the full set of one-particle matrix elements. 

We show an example of the Wigner function in Figs. 4(a) and 4(b) for a repulsive interaction. 
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Fig. 4: Friedel oscillations in phase space: Wigner function for N — 1000 interacting spinless fermions at zero 
temperature in the one-dimensional harmonic trap. In (a), the cross section W(z,k = 0) is plotted versus 
dimcnsionless distance z/l displaying the region around the classical turning point Lp = l^/2N — 1. 1 is the 
oscillator length. In (b), the corresponding function W(z — 0,k) is shown versus dimensionless momentum 
kl. Note the increase of the oscillation amplitude near the classical boundary. The plot refers to a repulsive 
interaction with K — 1/3. Repulsive interactions enhance Friedel oscillations in the spatial direction and 
suppress them in the momentum direction. 

The amplitudes of the Friedel oscillations increase near the classical turning point. This was discussed for the 
non-interacting case |41| . 

Neglecting non-diagonal matrix elements in Eq. (|65fl would give W(z, k = 0) = W(z = 0, k — > z/l 2 ). The significant 
differences between Fig. 4(a) and Fig. 4(b) thus stem from the non-diagonal matrix elements, which are particular 
relevant for repulsive interactions. 

The static pair correlation function with respect to the center of the trap 



C(z, z' = 0) = (ft(z) xj){z' = 0)) = 



2- 



dke-' kz W{z/2,k), 



(69) 



as well as any other one-particle property of the many-particle system are also contained in the Wigner function. 

The central pair correlation function l|69|) is displayed in Fig. 5. It is noted that the wave length of the intrinsic 
periodicity in the center of the trap is twice that of the Friedel oscillations, i.e., A = 2it /kp. In fact, the central pair 
correlation function for non-interacting fermions is given by 



C (z,z' = 0) = 



M 



n=0 



H 2n (z/l) 
4"n! 



sin(kpz) 



(70) 



provided N = 2M is large and z is restricted to the center of the trap (\z\ <C Lp). We expect that interactions modify 
Eq. (|70|l according to 



C(l < z < L F ,z' = 0) oc 



sin(fc^z) 
l{z/l) a c 



(71) 



i.e., the anomalous dimension of the correlation function is ac = 1 + 270 = (K + l/K)/2, as in a Luttinger liquid. 
This is confirmed numerically by displaying the corresponding envelope in Fig. 5. 
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Fig. 5: Dimensionless correlation functions Cq{z, z' = 0)1 and C(z,z' = 0)1 in units of the inverse of the oscillator 
length I versus dimensionless distance z/l from the center of the trap for N = 1000 spinless fermions at zero 
temperature. Thin curve shows non-interacting function Co while the thick curve gives C for an attractive 
interaction with K — 3. Oscillations have twice the period of Friedel oscillations. Envelope to the full curve 
displays power law decay with anomalous dimension (K + l/K)/2. 

The correlation function C of the interacting system thus decays faster than Co and this effect is invariant under 
K -> l/K. 

VI. SUMMARY 

The fact that backscattering dominates the interaction between identical one-dimensional fermions confined to a 
harmonic trap makes it necessary to include backscattering into the bosonization method. This can be done exactly 
by supplementing the squeezing transformation (|16f) with the displacement transformation (|18fl . This changes (in fact 
complicates) the result H40(l for the one-particle matrix elements by altering the argument in the sine-function. This 
modification destroys the symmetry 

{c\ N _ x _ n _ p C2N-l-n+p) = <5p,0 - (cl_ p C n+p ). (72) 

In order to evaluate the new result, we introduce in analogy to the Luttinger model a simplified form for the 
interaction coefficients which uses only one coupling constant K. K = 1 corresponds to the non-interacting case 
studied in 0, liBl l46l |47| . K > 1 corresponds to attraction. For values of K near unity, the Fermi edge at zero 
temperature is already significantly smoothed out by the interactions, an effect which becomes more pronounced for 
stronger coupling. 

We then study particle and momentum densities for various coupling strengths. Friedel oscillations in phase space are 
characteristic for these quantities, which extend progressively into the classically forbidden region when the coupling 
strength is increased. The results confirm an earlier finding that attractive interactions decrease the amplitude of 
the Friedel oscillations in real space while repulsion enhances it. The Friedel oscillations in momentum space behave 
oppositely. 

Finally, the central pair correlation function is calculated. The basic period of the pair correlation function in 
the center of the trap is 27r/fc^ in contrast to n/kp of the Friedel oscillations. Interactions do not affect these basic 
periodicities: Even the increase of the Friedel period in the particle density near the classical boundaries is correctly 
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given by the non-interacting theory. Interactions, however, modify the power law decay of the correlation function 
inside the trap. Our numerical results are consistent with the value 

K + l/K 

a c = (73) 

for the anomalous dimension ac of the correlation well inside the trap, in accordance with the predictions for a 
Luttinger liquid. 

We also assessed the validity of the bosonization method, comparing its result with those of direct numerical 
diagonalization in fermionic Hilbert space. 
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VII. APPENDIX 
A model for the interaction coefficients 

The result (|4(J|) for the one-particle matrix elements contains three sets of interactions coefficients: {a m }, {7 m }, 
and {V c (m)}. The first two sets depend on the basic interaction coefficients V c (m) via the parameter £ TO according to 
Eq. (|20|l . The interaction coefficients V c (m) appear also directly in W%. In order to get explicit results for the matrix 
elements, one must be able to perform the summations in W\ and W^. This requires models for the m-dependence of 
the above interaction coefficients. 

Defining V c = V c /htO£, an explicit form for a m is 

am = / = -T7? • ( A1 ) 

201 - V c (m)] 2 4A ™ 

Similarly, the central coupling constants K m determine "f m and the renormalized oscillator energies e m according 

(1 - K m f 2K m 
lm= AK m ' em = nWt KlT^- (A2) 

Following the procedure in the Luttinger model, we adopt exponential decays, thus we make the ansatz: 

a m = a exp(-r Q m/2), (A3) 

and 

7m = 7oexp(-r 7 m). (A4) 
Since the signs of V c (m) do not depend on m, a m is related to 7 m by 



to 



a m = -sgn(V c ) V7m(i + j m ), (A5) 

and one finds 



-sgn(V r c ) exp(-r 7 m/2)y / 7o[l + 7o exp(-r 7 m)]. (A6) 
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Assuming that r 7 is very small, < r 7 < 1, and that the relevant values of m obey m < 1/V 7 , Eq. (|A6|I leads to 

a m exp(-r 7 m/2) a , (A7) 



with 



a = -sgn(F c ) VToCl + To). (A8) 

Thus, we can set 

r Q = r 7 = r. (A9) 

Another set of coupling parameters involves directly the backscattering coefficients V c (m) via in Using 
Eqs. (EH) and (J2SJ), Eq. becomes 



V c (m) 



4m [1 - V c (2m)] 

We take advantage of the slow decay of interaction coefficients with m and write 



with another decay constant r c <C 1 and with 



(A10) 



6m = —j—Vcejf exp(-r c m), (All) 
4m 



" l-K(l) 2 
This gives a useful result for the backscattering function (|^Tt 

C(u) = -Y^H V" 1 exp[(-m(2ro + r c )l = m h _ cxp (_ rc _ 2iu)]. (A13) 

4 ^— ' m "4 

m— 1 

It is noted that consistency requires \V C \ < 1 and K > 0. 

Comparing a(m) in Eq. (|A1|) with Eqs. (|A10I) and IjAlljl . leads to 

a m = -V ceff exp(-r c m) { 1 U c (2m) = . 



2^1 -k 2 M 

Suppressing the weak m-dependence in the curly brackets, i.e., 

V c {m) -> (A15) 

allows the identifications 

r.-^-jj- (AW) 

Furthermore, using 
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gives 



The main coupling constant 



ao 



Vc 

2Jl - V 2 ' 



V< 



ceff 



(A17) 




(A18) 



K 



l + V c 



, i-v c 

determines all relevant couplings and the renormalized energy according to 



(A19) 



Vc = 



K 2 



1 



K 2 



1 



K 2 + l 



a = 



7o = 



(K-lf 



AK 



= TlLOf 



2K 



K 2 + l 



(A20) 



Only the two parameters K and r = r a = r~ i = 2r c remain. K is determined by the physical backscattering 
strength V c (l) = V c according to Eq. JXHJ). In {2|j, it was argued that r = 1/ y/N is a reasonable choice for the 
decay constant. Finally, it is noted that interactions always decrease the renormalized excitation energy e below the 
non-interacting value hu>£. 
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